##########################################################################################

library(data.table)
library(optparse)
library(WGCNA)
library(parallel)
library(dplyr)
library(RColorBrewer)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--MEDissThres"), type = "character"),
    make_option(c("--minClusterSize"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/wgcnv/CombineCounts.FilterLowExpression-MergeMutiSample.varianceStabilizingTransformation.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/wgcnv"
    # Merge close modules
    MEDissThres=0.25
    minClusterSize=100

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
MEDissThres <- as.numeric(opt$MEDissThres)
minClusterSize <- as.numeric(opt$minClusterSize)

dir.create(out_path , recursive = T)

##########################################################################################
# The following setting is important, do not omit.
options(stringsAsFactors = FALSE);
# Allow multi-threading within WGCNA. This helps speed up certain calculations.
# At present this call is necessary.
# Any error here may be ignored but you may want to update WGCNA if you see one.
# Caution: skip this line if you run RStudio or other third-party R environments.
# See note above.
enableWGCNAThreads(nThreads=10)
# Load the data saved in the first part

## I am not aware of a principle from which one could derive an “appropriate” value, but in practice, 
## on data sets with 50-100 samples, using 0.15 to 0.2 has worked fairly well. 
## For fewer samples a larger valuse (0.25 to 0.3) may be warranted. 
## If you want larger modules, increase the value; 
## if you want smaller modules at the risk of having redundant modules 
## (modules with very similar functional annotation and very similar fuzzy module membership), 
## you can decrease the value to say 0.10, maybe even lower if you have lots of samples (hundreds or more).

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))

##########################################################################################
## https://mp.weixin.qq.com/s/2VXprzgJSYJ39AhZ_FprCQ
## 1至少20个样本以上，越多越好；
## 2可以过滤点低表达或者低方差的基因，以减少干扰信息。但不太建议直接使用差异基因。
## 3WGCNA最初用于芯片测序数据，也适用于RNA-seq数据。关于RNAseq标准化，由于不涉及到不同基因之间的比较，所有常规标准化方式都可以。
## https://www.homedt.net/233619.html
## 输入数据形式如果有批次效应,需要先进行去除；
## 处理RNAseq数据，需要采用DESeq2的varianceStabilizingTransformation方法，或将基因标准化后的数据（如FPKM、CPM等）进行log2(x+1)转化
## https://www.biostars.org/p/280650/

##########################################################################################
##  选择合适的软阈值β
exp_dat <- dat_tpm[,1:(ncol(dat_tpm)-1)]
rownames(exp_dat) <- dat_tpm$gene_id
exp_dat <- t(exp_dat)

powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(exp_dat, powerVector = powers, verbose = 5)

out_name <- paste0( out_path , "/sft_plot.pdf" )

pdf(out_name)
cex1 = 0.9
par(mfrow = c(2,1))
### （1）是否符合幂律分布；
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red")
abline(h=0.80,col="red")
### （2）节点的平均连接度
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
#par(mfrow = c(1,1))
dev.off()

## 基本在6时，拟合幂律分布的结果是比较好的；同时节点的凭据连通性也趋于稳定了。之后会使用这个参数建立网络
## To maximize strong correlations between genes, we assigned a soft-thresholding power β = 8 since it was the 
## lowest β with the highest R2 value that passed the 0.8 threshold 
power <- sft$fitIndices[,1][-sign(sft$fitIndices[,3])*sft$fitIndices[,2] > 0.8][1]
#power <- 6

##########################################################################################
## 挑选模块Hub基因应该考虑
## (1) 模块内基因连接度
# WGCNA offers a function "intramodularConnectivity" that actually also computes the total connectivity (within the whole network). 
# In addition to total and intramodular connectivity it also gives the extra-modular connectivity (kTotal-kWithin), 
# and the difference of the intra- and extra-modular connectivities (kIn-kOut) for all genes. 
# The value returned from the function is a data frame with one column for each of the measures mentioned above. 
# Intramodular connectivity is calculated by summing adjacency entries (excluding the diagonal) to other nodes within the same module. 
# Intramodular connectivity (kWithin)
# After raising the module membership to a power of 6, it is highly correlated with the intramodular connectivity (kWithin).
## https://www.jianshu.com/p/c9c73e1093fd
## 初步构建gene module
## 以TOM dissimilarity矩阵作为输入，进行聚类分析。
##########################################################################################

adjacency = adjacency(exp_dat, power = power)

# Turn data expression into topological overlap matrix
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM

# Plot gene tree
geneTree = hclust(as.dist(dissTOM), method = "average");   # 用于后续cutreeDynamic()，对gene tree进行裁剪

out_name <- paste0( out_path , "/Gene_cluster.pdf" )
pdf(file = out_name, width = 12, height = 9);
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04);
dev.off()

##########################################################################################
## 使用cutreeDynamic()进行聚类分析的优化
## Module identification using dynamic tree cut
## https://jneuroinflammation.biomedcentral.com/articles/10.1186/s12974-019-1433-4#Sec2
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = FALSE , method = "tree" , minClusterSize = minClusterSize )

# Convert numeric labels into colors
dynamicColors = labels2colors(dynamicMods)
#table(dynamicColors)
# Plot the dendrogram and colors underneath
out_name <- paste0( out_path , "/Module_tree.pdf" )
pdf(file = out_name, width = 8, height = 6);
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",dendroLabels = FALSE,
                    hang = 0.03,addGuide = TRUE, guideHang = 0.05,main = "Gene dendrogram and module colors")
dev.off()


merge = mergeCloseModules(exp_dat, dynamicColors, cutHeight = MEDissThres, verbose = 3) 
mergedColors = merge$colors  
mergedMEs = merge$newMEs  
# Plot merged module tree
out_name <- paste0( out_path , "/Merged_Module_Tree.pdf" )

pdf(file = out_name , width = 12, height = 9)  
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), 
                    c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = FALSE, 
                    hang = 0.03, addGuide = TRUE, guideHang = 0.05)  
abline(h=MEDissThres, col = "red")
dev.off()

##########################################################################################
## 将具体描述模块或基因与表型的关联
group <- sapply( strsplit(rownames(exp_dat) , "_" ) , "[" , 2 )
group <- factor(group , levels = unique(group)) #根据样本信息进行分组。

datTraits <- data.frame(samples=rownames(exp_dat),subtype=group)

design <- model.matrix(~0 + datTraits$subtype)
colnames(design) <- levels(factor(datTraits$subtype))

#计算合并后的模块的第一主成分ME。
MEs0 <- moduleEigengenes(exp_dat, mergedColors)$eigengenes
MEs <- orderMEs(MEs0); #不同颜色的模块的ME值矩阵。
moduleTraitCor = cor(MEs, design , use = "p") #计算不同模块和样本性状的相关性。
nSamples <- nrow(exp_dat)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples) #计算不同模块和样本性状相关性的p值。

# 可视化module与性状的相关性和它们的p值
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "");

dim(textMatrix) = dim(moduleTraitCor)

out_name <- paste0( out_path , "/group_Modules_heatmap.pdf" )
pdf(file = out_name , width = 15, height = 13)  
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(design),
xLabelsAngle = 0,
yLabels = row.names(moduleTraitCor),
ySymbols = row.names(moduleTraitCor),
colorLabels = TRUE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 1,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
dev.off()


##########################################################################################
softpower <- power
datExpr <- exp_dat

# Define numbers of genes and samples
nGenes = ncol(exp_dat);
nSamples = nrow(exp_dat);

## 1、计算module membership
## 使用的是WGCNA自带cor函数，使用皮尔逊计算相关性
MEs0 = moduleEigengenes(exp_dat, mergedColors)$eigengenes
MEs = orderMEs(MEs0)
geneModuleMembership = as.data.frame(cor(exp_dat, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
modNames = substring(names(MEs), 3) # 去除ME前缀
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");

## 2、计算Gene significance GS 的计算
geneTraitSignificance = as.data.frame(cor(datExpr, design, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.",names(geneTraitSignificance), sep="");
names(GSPvalue) = paste("p.GS.",names(geneTraitSignificance), sep="");

## 3、 计算连接度
## 跟TOM一样
connectivity <- abs(cor(datExpr,use="p"))^softpower
## 计算每个基因module内的连接度和module外的连接度以及total连接度。
Alldegrees <- intramodularConnectivity(connectivity, mergedColors) 
datKME <- signedKME(datExpr, MEs, outputColumnName="kME_MM.")#计算MM。
GS <- cor(datExpr,design,use="p") #注意，如果是design是离散变量，一定要将design放在后面

## 合并计算的自由度和KME
combin_data <- cbind(Alldegrees , datKME , GS)
combin_data$class <- mergedColors


## 计算每个module表达最相似的5%基因定义为hub
result_dat <- c()
for( type in unique(mergedColors) ){
     use_data <- subset( combin_data , class==type)
    
     ## membership找hub基因
     t_th <- as.numeric(quantile(use_data[,paste0("kME_MM." , type)] , seq(0,1,0.05))[20])
     use_data$Hub_membership <- ifelse(use_data[,paste0("kME_MM." , type)] >= t_th , "TRUE" , "FALSE")
     use_data <- use_data[order( use_data[,paste0("kME_MM." , type)] , decreasing=T),]
     use_data$Hub_membership_rank <- order(use_data[,paste0("kME_MM." , type)] , decreasing=T)/nrow(use_data)

     ## kwithin找hub基因
     t_th <- as.numeric(quantile(use_data$kWithin , seq(0,1,0.05))[20])
     use_data$Hub_kWithin <- ifelse(use_data$kWithin >= t_th , "TRUE" , "FALSE")
     use_data <- use_data[order(use_data$kWithin , decreasing=T),]
     use_data$kWithin_rank <- order(use_data$kWithin , decreasing=T)/nrow(use_data)

     result_dat <- rbind( result_dat , use_data )
}

result_dat$gene_id <- rownames(result_dat)

out_name <- paste0(out_path , "/wgcnv_HubModule.tsv")
write.table( result_dat , out_name , row.names = F , sep = "\t" , quote = F )

#### https://github.com/halryd/high_S100A_hub_genes/blob/master/08_chap3_select_hub_genes.Rmd
## In the tutorials of WGCNA they indicate three different way to measure the degree to which a gene is a hub. 
## It is shownd that GS and MM correlate quit well an that. MM and 
## 1. Gene significance (GS)
## 2. Module membership (MM)
## 3. Module connectivity (MC?)


## - In paper "Inflammatory, regulatory, and autophagy co-expression modules and hub genes underlie the peripheral immune response 
## to human intracerebral hemorrhage" they are using top 5% highest membership to their respective module as definition of a hub gene-
## References for 5% membership: @RN7


## In "Gene co-expression network analysis reveals common system-level properties of prognostic genes across cancer types" they  defined the 1% (or 5%) of nodes with the highest connectivity as hubs 25–27.
## They write further "Given a network, we then obtained several key network properties such as the edge weight, node connectivity and modularity. Connectivity was defined as the sum of the weights across all the edges of a node, and the top 1% (or 5%) of the genes with the highest connectivity in the network were defined as hub genes."
## References for 5% connectivity @RN5

